Tarea 3

Componentes principales

Pregunta 1

Si dos variables $X$ y $Y$ tienen covarianza $S = \begin{pmatrix} a && b \\ c && d \end{pmatrix}$, muestre que cuando $c \neq 0$, la primera componente principal está dada por \[ \sqrt{\frac{c^2}{c^2+(V_1-a)^2}}X+\frac{c}{|c|}\sqrt{\frac{(V_1-a)^2}{c^2+(V_1-a)^2}}Y \] donde $V_1$ es la varianza explicada por la primera componente principal.

Pregunta 3

Considere los datos en el archivo T8-5.DAT correspondientes a un tramo censal. Suponga que las observaciones en la variable $X5$ = valor de la mediana de hogares fue registrada en unidades de diez miles más que de cientos de miles de dólares (es decir, multiplique todos los datos listados en la sexta columna por 10).

En cada caso, compare las estimaciones con los datos en diez miles y cientos de miles (son dos matrices de covarianzas) para las componentes pincipales, trate de obtener una interpretación y de explicar el efecto de cambiar de escala.

col_names_censo <- c('poblacion', 'p_profesional', 'p_empleados_m16',
                     'p_empleados_gob', 'mediana_hogares')
censo <- read_table('./data/T85.DAT', col_names=col_names_censo)
rm(col_names_censo)

censo_escala2 <- censo %>%
  mutate(mediana_hogares = 10*mediana_hogares)

pc_orig <- prcomp(censo, scale=F) %>%
  print()

## Standard deviations (1, .., p=5):
## [1] 2.6326932 1.3360929 0.6242194 0.4790918 0.1189747
## 
## Rotation (n x k) = (5 x 5):
##                         PC1         PC2          PC3         PC4
## poblacion       -0.78120807  0.07087183 -0.003656607  0.54171007
## p_profesional   -0.30564856  0.76387277  0.161817438 -0.54479937
## p_empleados_m16 -0.33444840 -0.08290788 -0.014841008  0.05101636
## p_empleados_gob -0.42600795 -0.57945799 -0.220453468 -0.63601254
## mediana_hogares  0.05435431  0.26235528 -0.961759720  0.05127599
##                          PC5
## poblacion        0.302039670
## p_profesional    0.009279632
## p_empleados_m16 -0.937255367
## p_empleados_gob  0.172145212
## mediana_hogares -0.024583093

summary(pc_orig)

## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5
## Standard deviation     2.6327 1.3361 0.62422 0.47909 0.11897
## Proportion of Variance 0.7413 0.1909 0.04168 0.02455 0.00151
## Cumulative Proportion  0.7413 0.9323 0.97394 0.99849 1.00000

pc_nuevo <- prcomp(censo_escala2, scale=F) %>%
  print()

## Standard deviations (1, .., p=5):
## [1] 7.1392827 2.5787380 1.1916363 0.4793601 0.1190091
## 
## Rotation (n x k) = (5 x 5):
##                         PC1         PC2         PC3          PC4
## poblacion        0.05765977 -0.78220765 -0.02307185  0.541282426
## p_profesional   -0.03281610 -0.35032828 -0.76422159 -0.540439691
## p_empleados_m16  0.03467366 -0.32717525  0.10109097  0.051177198
## p_empleados_gob  0.07557524 -0.39085238  0.63207658 -0.642117054
## mediana_hogares -0.99432619 -0.07491367  0.07545108  0.002204195
##                          PC5
## poblacion       -0.302171422
## p_profesional   -0.009137870
## p_empleados_m16  0.937504990
## p_empleados_gob -0.172301138
## mediana_hogares  0.002375237

summary(pc_nuevo)

## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5
## Standard deviation     7.1393 2.5787 1.19164 0.47936 0.11901
## Proportion of Variance 0.8598 0.1122 0.02395 0.00388 0.00024
## Cumulative Proportion  0.8598 0.9719 0.99589 0.99976 1.00000

En la escala original, la primera componente principal explica el 74% de la varianza. Le asigna el coeficiente de más peso a poblacion, pero salvo por mediana_hogares, las otras variables tienen pesos moderados y similares.

Al hacer el cambio de escala en mediana_hogares, sin embargo, la primera componente se vuelve casi mediana_hogares por si sola, pues los demás coeficientes son casi cero. Esto puede explicarse porque componentes principales está pensado para explicar la mayor varianza posible, y en esta nueva escala, como $\textrm{Var}(10X_6)=100\textrm{Var}(X_6)$, la varianza aumenta significativamente, y también lo hace en comparación con las demás variables, como podemos ver usando

apply(as.matrix(censo), 2, sd)

##       poblacion   p_profesional p_empleados_m16 p_empleados_gob 
##       2.0754652       1.3294632       0.8948008       1.4033798 
## mediana_hogares 
##       0.7101973

apply(as.matrix(censo_escala2), 2, sd)

##       poblacion   p_profesional p_empleados_m16 p_empleados_gob 
##       2.0754652       1.3294632       0.8948008       1.4033798 
## mediana_hogares 
##       7.1019731

Pregunta 5

Consideren la matriz de correlaciones dada abajo. Los datos originales corresponden a las mediciones de 8 variables de química sanguínea de 72 pacientes en un estudio clínico. (Jolliffe, 2002). La matriz de correlaciones de las variables rblood, plate, wblood, neut, lymph, bilir, sodium y potass, en ese orden, está dada en $S$, y la desviación estándar de cada variable en $\mathbf{\sigma}$.

Aplique componentes principales a la matriz de covarianzas y a la matriz de correlaciones y explique las diferencias. Basado en la observación anterior, ¿sobre qué debería hacerse el análisis?

R <- read_csv('./data/5-S.csv')
sigma <- c(0.371, 41.253,  1.935,  0.077,  0.071,  4.037,  2.732,  0.297)
names(sigma) <- names(R)

S <- t((t(as.matrix(R) * sigma)*sigma))

pc_S <- princomp(covmat=as.matrix(S), cor=F) 
pc_R <- princomp(covmat=as.matrix(R), cor=T)

pc_S$loadings

## 
## Loadings:
##        Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
## rblood                              0.943  0.329              
## plate  -0.999                                                 
## wblood        -0.192        -0.981                            
## neut                                              0.758  0.650
## lymph                                            -0.649  0.760
## bilir          0.961  0.195 -0.191                            
## sodium         0.193 -0.979                                   
## potass                              0.329 -0.942              
## 
##                Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
## SS loadings     1.000  1.000  1.000  1.000  1.000  1.000  1.000  1.000
## Proportion Var  0.125  0.125  0.125  0.125  0.125  0.125  0.125  0.125
## Cumulative Var  0.125  0.250  0.375  0.500  0.625  0.750  0.875  1.000

summary(pc_S)

## Importance of components:
##                            Comp.1      Comp.2     Comp.3      Comp.4
## Standard deviation     41.2877486 3.880212624 2.64197339 1.624583979
## Proportion of Variance  0.9856182 0.008705172 0.00403574 0.001525986
## Cumulative Proportion   0.9856182 0.994323381 0.99835912 0.999885108
##                             Comp.5       Comp.6       Comp.7       Comp.8
## Standard deviation     0.353951757 2.561722e-01 8.510631e-02 2.372715e-02
## Proportion of Variance 0.000072436 3.794288e-05 4.187837e-06 3.255049e-07
## Cumulative Proportion  0.999957544 9.999955e-01 9.999997e-01 1.000000e+00

pc_R$loadings

## 
## Loadings:
##        Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
## rblood -0.194  0.417  0.400  0.652  0.175 -0.363  0.176  0.102
## plate  -0.400  0.154  0.168        -0.848  0.230 -0.110       
## wblood -0.459         0.168 -0.274  0.251  0.403  0.677       
## neut   -0.430 -0.472 -0.171  0.169  0.118        -0.237  0.678
## lymph   0.494  0.360        -0.180 -0.139  0.136  0.157  0.724
## bilir   0.319 -0.320 -0.277  0.633 -0.162  0.384  0.377       
## sodium  0.177 -0.535  0.410 -0.163 -0.299 -0.513  0.367       
## potass  0.171 -0.245  0.709         0.198  0.469 -0.376       
## 
##                Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
## SS loadings     1.000  1.000  1.000  1.000  1.000  1.000  1.000  1.000
## Proportion Var  0.125  0.125  0.125  0.125  0.125  0.125  0.125  0.125
## Cumulative Var  0.125  0.250  0.375  0.500  0.625  0.750  0.875  1.000

summary(pc_R)

## Importance of components:
##                           Comp.1    Comp.2    Comp.3     Comp.4     Comp.5
## Standard deviation     1.6710100 1.2375848 1.1177138 0.88227419 0.78839505
## Proportion of Variance 0.3490343 0.1914520 0.1561605 0.09730097 0.07769584
## Cumulative Proportion  0.3490343 0.5404863 0.6966468 0.79394778 0.87164363
##                            Comp.6     Comp.7     Comp.8
## Standard deviation     0.69917350 0.66002394 0.31996216
## Proportion of Variance 0.06110545 0.05445395 0.01279697
## Cumulative Proportion  0.93274908 0.98720303 1.00000000

Las componentes principales, por construcción, están hechas para explicar el mayor porcentaje de la varianza. Por lo tanto, cuando se utiliza $S$, como una variable (plate) tiene muchísima más variabilidad que las demás, la primera componente principal le asignará un coeficiente muy grande.

Aún más notorio es al comparar la proporción de la varianza que explica cada componente: cuando se utilizan las covarianzas, la primer componente principal explica casi toda. Peor: para explicar la misma proporción con el modelo de correlaciones se necesitan siete de las ocho componentes.

El análisis usando $S$ sólo dice lo que podíamos ver desde $\mathbf{\sigma}$: que plate varía mucho, y en este sentido es mejorel de $R$. Sin embargo, hay que ser cuidadosos al interpretar porque los de $R$ son de las variables estandarizadas, que están en desviaciones desde la media y no las medidas originales.

Regresión lineal

Pregunta 7

Considere los datos sobre venta de casas en T7-1.DAT. Este archivo tiene las siguientes variables:

$X_1$ = tamaño de construcción (en 100 $\textrm{ft}^2$)
$X_2$ = precio de valuación (en miles de USD)
$Y$ = precio de venta (en miles de USD)

Estime el modelo de regresión lineal $\mathbf{Y} = \mathbf{\beta'x}+\mathbf{\varepsilon}$. Haga el diagnóstico del modelo sobre los residuales.

casas <- read_csv('./data/t71.csv')
casas_lm <-lm(p_venta~tamaño+p_valuacion, data=casas) 

casas_tidy <- tidy(casas_lm)
casas_glance <- glance(casas_lm)

Parece ser que la valuación anterior no es estadísticamente significativa. Además, $R^2=0.8344$ y $R^2_{\textrm{aj}}=0.8149$, por lo que el modelo explica buena parte de la varianza.

Primero queremos buscar no-linealidad en los datos a través de una gráfica de residuales contra valores ajustados:

studentized <- rstudent(casas_lm)
casas_lm <- augment(casas_lm) %>%
  mutate(.student.resid = studentized) 
rm(studentized)

casas_lm %>%
  ggplot(aes(x=.fitted, y=.resid)) +
  geom_point() +
  geom_smooth() +
  geom_smooth() +
  labs(x='Fitted values', y='Residuals')

casas_lm %>%
  ggplot(aes(x=.fitted, y=sqrt(abs(.std.resid)))) +
  geom_point() +
  geom_smooth() +
  ggtitle('Scale-location') +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x='Fitted values', y=expression(sqrt(.std.resid)))

Aunque la LOESS parece tener forma de u en cerca de 75,pero en general la dispersión de los puntos no parece seguir patrones y la varianza se ve más o menos constante.

Queremos concentrarnos en un valor que puede ser problemático porque se ve lejos de los demás. Aprovechando que son sólo dos variables, podemos hacer una gráfica de dispersión para corroborar que efectivamente está lejos de los demás puntos.

casas_lm %>%
  ggplot(aes(x=tamaño, y=p_valuacion)) +
  geom_point()+ 
  ggtitle('Predictor space') +
  theme(plot.title = element_text(hjust = 0.5))

Para medir la influencia del punto aunque sea visualmente,

casas_lm %>%
  ggplot(aes(x=1:20, y=.cooksd)) +
  geom_point() +
  geom_line() +
  labs(x='Observation number', y="Cook's Distance")

casas_lm %>%
  ggplot(aes(x=1:20, y=.hat)) +
  geom_point() +
  geom_line()+
  labs(x='Observation number', y='Leverage')

p <- casas_lm %>%
  ggplot(aes(x=.hat, y=.student.resid)) +
  geom_point(aes(size=.cooksd)) +
  labs(title='Influence plot', 
       x='Leverage', 
       y='Studentized residuals', 
       size="Cook's Distance") +
  theme(plot.title=element_text(hjust=0.5))

ggplotly(p)

La distancia de Cook indica que la observación 16 podría tener mucha influencia sobre el modelo, y el apalancamiento confirma que es por estar lejos de las demás en el espacio de variables. Sin embargo, en la gráfica de influencia no se ve que el el valor de la respuesta esté muy lejos de la tendencia.

Para buscar correlación entre los errores,

casas_lm %>%
  ggplot(aes(x=1:20, y=.resid)) +
  geom_line() +
  geom_point() +
  labs(x='Observation number', y='Residuals')

Acf(casas_lm$.resid)

Y parece que no hay correlación significativa entre errores. Para verificar que se distribuyan normal y poder hacer inferencia tranquilamente, hacemos una qqplot

casas_lm %>%
  ggplot(aes(sample=.student.resid)) +
  geom_qq_band(alpha=0.5, distribution='norm') +
    #geom_qq_band(alpha = 0.5, distribution='t', dparams=17) +
  stat_qq_line() +
  stat_qq_point() +
  labs(x = "Theoretical quantiles", y = "studentized residuals")

Aunque los puntos no forman una línea muy recta en las colas, sólo uno se sale de los intervalos de confianza, lo cual es normal para una prueba al 5% de significancia sobre 20 datos. Nos damos por satisfechos suponiendo normalidad.

Obtenga un intervalo de predicción del 95% para el precio de venta $Y_0$ usando $X_1 =17$ y $X_2 =46$.

x0 <- c(1, 17, 46)
X <- casas %>%
  mutate(const=1) %>%
  select(const, tamaño, p_valuacion) %>%
  as.matrix()
Y <- casas$p_venta

s2 <- casas_glance$sigma^2
beta_gorro <- casas_tidy$estimate

exts <- qt(0.95/2, 20-2-1, lower.tail = FALSE)*sqrt(s2*(1+t(x0)%*%solve(t(X)%*%X)%*%x0))

x0%*%beta_gorro-exts

##          [,1]
## [1,] 77.41626

x0%*%beta_gorro+exts

##          [,1]
## [1,] 78.24337

Realice una prueba LRT sobre $ H_0 : _2 = 0$ con nivel de significancia de $\alpha= 0.05$

modelo_H0 <- lm(p_venta ~ tamaño, data=casas)
modelo_Ha <- lm(p_venta ~ tamaño+p_valuacion, data=casas)

anova(modelo_H0, modelo_Ha)

## Analysis of Variance Table
## 
## Model 1: p_venta ~ tamaño
## Model 2: p_venta ~ tamaño + p_valuacion
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1     18 205.3                           
## 2     17 205.0  1    0.3027 0.0251  0.876

glance(modelo_Ha)

## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
## *     <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>  <dbl> <dbl> <dbl>
## 1     0.834         0.815  3.47      42.8 2.30e-7     3  -51.7  111.  115.
## # ... with 2 more variables: deviance <dbl>, df.residual <int>

F_stat <- ((205.3-205)/(18-17))/glance(modelo_Ha)$sigma^2 #(Leí los numeritos de las llamadas de arriba)
p_value <- 1-pf(F_stat, 1, 17)
p_value

## [1] 0.8765281

No hay evidencia para rechazar la hipótesis nula, por lo que nos quedamos con el modelo que no incluye la valuación pasada.

Obtenga la gráfica de la $C_p$ de Mallows para este problema.

cp <- regsubsets(p_venta~tamaño+p_valuacion, data=casas) %>%
  summary() %>%
  '$'(cp)
data_frame(p=1:2, cp=cp) %>%
  ggplot(aes(p,cp)) +
  geom_point() +
  geom_abline(aes(slope=1,intercept=0))

Una vez más, indica que nos quedemos con el modelo más pequeño.

Pregunta 9

El Berkeley Guidance Study fue un estudio hecho a través del tiempo para seguir a un grupo de niños y niñas desde que nacieron en Berkeley entre enero de 1928 y junio de 1929, hasta al menos los 18 años. Los datos se encuentran en el archivo BGSall.DAT.

Para las niñas, obtengn las estadísticas sumarias usuales (medias, desviaciones estándar y correlaciones) para todas las variables, excepto Case y Sex. Obtenga la matriz de scatterplots para las variables de la edad 2, las variables de edad 18 y HT18. Resuma la información que se pueda extraer de la gráfica.

colnames <- c('Sex', 'WT2', 'HT2', 'WT9', 'HT9', 'LG9', 'ST9', 
              'WT18', 'HT18', 'LG18', 'ST18', 'soma', 'case')
berkeley <- read.table('./data/BGSall.DAT', skip=24, col.names=colnames) %>%
  as_data_frame()
rm(colnames)

women_mean <- berkeley %>% 
  filter(Sex==1) %>%
  select(-Sex, -soma) %>%
  colMeans()

women_cor <- berkeley %>% 
  filter(Sex==1) %>%
  select(-Sex, -soma) %>%
  cor()

women_sd <- berkeley %>% 
  filter(Sex==1) %>%
  select(-Sex, -soma) %>%
  as.matrix() %>%
  apply(2, sd)

ggpairs(data = berkeley,
        mapping=ggplot2::aes(colour=as.factor(Sex)),
        columns=2:12,
        title='Berkeley Study Group'
        )

Hay correlación positiva entre HT18 y casi todas las variables, particularmente la estatura a los dos años y a los nueve.

Ajuste el modelo de regresión lineal \[ \textrm{HT18}|\mathbf{X} = \alpha_0 + \alpha_1\textrm{WT2} + \alpha_2\textrm{HT2} + \alpha_3\textrm{WT9} + \alpha_4\textrm{HT9} + \alpha_5\textrm{LG9} + \alpha_6\textrm{ST9} + \epsilon \] suponiendo que la varianza es constante y dé los estimados y los errores estándar para todos los parámetros, así como el valor del coeficiente de determinación.

bk_model <- lm(HT18~WT2+HT2+WT9+HT9+LG9+ST9, data=berkeley)
tidy(bk_model)

## # A tibble: 7 x 5
##   term        estimate std.error statistic     p.value
##   <chr>          <dbl>     <dbl>     <dbl>       <dbl>
## 1 (Intercept)  41.3      21.7        1.91  0.0590     
## 2 WT2           1.42      0.494      2.87  0.00478    
## 3 HT2           0.127     0.255      0.496 0.621      
## 4 WT9          -0.0898    0.270     -0.333 0.740      
## 5 HT9           0.969     0.177      5.46  0.000000236
## 6 LG9          -1.20      0.569     -2.12  0.0360     
## 7 ST9           0.0994    0.0433     2.29  0.0235

glance(bk_model)$adj.r.squared

## [1] 0.5124361

Muestre numéricamente que $R^2$ es el mismo valor que el coeficiente de correlación entre HT18 y los valores ajustados de la ecuación de la pregunta anterior. Haga una gráfica de la respuesta contra los valores ajustados, y dé una interpretación visual del coeficiente de determinación.

r_squared <- glance(bk_model)$r.squared
bk_model_aug <- bk_model %>%
  augment()
corr_calc <- cor(bk_model_aug$HT18, bk_model_aug$.fitted)^2

r_squared/corr_calc

## [1] 1

bk_model_aug %>%
  ggplot(aes(HT18, .fitted)) +
  geom_point() +
  geom_smooth(method='lm', se = F)

El valor de $R^2$ dice que la recta explica cerca del 53% de la variabilidad en las observaciones, cosa que vemos en que los puntos están muy dispersos al rededor de la línea.

De las unidades de medición de cada uno de los coeficientes de regresión estimados.

$\alpha_1$, $\alpha_3$ y $\alpha_6$ no tienen unidades
$alpha_2$, $\alpha_4$ y $\alpha_5$ están en $\textrm{kgcm}^{-1}$

Obtenga pruebas de que cada coeficiente es igual a 0, y dar el valor de la estadística de prueba, su p−valor, y un breve resumen del resultado.

Haciendo las pruebas margianles:

summary(bk_model)

## 
## Call:
## lm(formula = HT18 ~ WT2 + HT2 + WT9 + HT9 + LG9 + ST9, data = berkeley)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.4766  -4.5606   0.1421   4.0278  11.7021 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 41.29577   21.67489   1.905  0.05898 .  
## WT2          1.41695    0.49351   2.871  0.00478 ** 
## HT2          0.12656    0.25501   0.496  0.62054    
## WT9         -0.08979    0.26981  -0.333  0.73984    
## HT9          0.96851    0.17740   5.460 2.36e-07 ***
## LG9         -1.20494    0.56873  -2.119  0.03604 *  
## ST9          0.09936    0.04333   2.293  0.02346 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.176 on 129 degrees of freedom
## Multiple R-squared:  0.5341, Adjusted R-squared:  0.5124 
## F-statistic: 24.65 on 6 and 129 DF,  p-value: < 2.2e-16

Los valores p indican que sólo podemos rechazar la hipótesis de coeficientes nulos con confianza superior a 95% en WT2, HT9, LG9 y ST9.

Probar la hipótesis de que los predictores a edad 2 no son necesarios en el modelo (i.e. prueba que $NH : \alpha_1 = \alpha_2 = 0$ versus la alternativa general).

anova(lm(HT18~WT9+HT9+LG9+ST9, data=berkeley),
      bk_model)

## Analysis of Variance Table
## 
## Model 1: HT18 ~ WT9 + HT9 + LG9 + ST9
## Model 2: HT18 ~ WT2 + HT2 + WT9 + HT9 + LG9 + ST9
##   Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
## 1    131 5380.8                                
## 2    129 4919.8  2    460.99 6.0437 0.003098 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La estadística F da suficiente evidencia para rechazar la hipótesis nula con 97% de confianza, indicando que es mejor el modelo con las variables a los dos años.

De intervalos de confianza del 95 % para $\alpha_4$ y para $\alpha_4 - \alpha_2$ de la ecuación.

Suponiendo normalidad y homocedasticidad, los coeficientes siguen \[ \mathbf{\hat{\alpha}} \sim \mathcal{N}(\mathbf{\alpha}, \sigma^2(X^\top X)^{-1}) \]

Por lo que podemos dar intervalos simultáneos

alpha <- bk_model$coefficients
s_alpha <- vcov(bk_model)

print('Intervalo para alpha 4')

## [1] "Intervalo para alpha 4"

lim_int <- sqrt(s_alpha[5,5])*sqrt((6+1)*pf(6+1,136-(6+1),0.05, lower.tail=F))
alpha[[4]] + c(-1,1)*lim_int

## [1] -0.5302151  0.3506350

print('')

## [1] ""

print('Intervalo para alpha 4 - alpha 2')

## [1] "Intervalo para alpha 4 - alpha 2"

a <- c(0, 0, -1, 0, 1, 0, 0)
centro <- a%*%alpha
lim_int <- sqrt(t(a)%*%s_alpha%*%a)* sqrt(((6+1)*(136-1)/(136-(6+1)))*pf(6+1,136-(6+1),0.05, lower.tail=F))
c(centro -lim_int, centro+lim_int)

## [1] -0.1317404  1.8156470

Pregunta 10

Con datos de alguna fuente de datos abiertos, proponer (a) un modelo de regresión múltiple, y (b) un modelo de componentes principales. En cada caso, realizar el análisis correspondiente, y elaborar un reporte del análisis en cada caso, de no más de una cuartilla en cada caso.

rm(list=ls())
power_plant <- read_excel('./data/CCPP/Folds5x2_pp.xlsx')
names(power_plant) <- str_replace_all(names(power_plant), ' ', '_')

pp_models <- regsubsets(PE~., data=power_plant, nvmax = 2500)

data_frame(p=1:4, cp=summary(pp_models)$cp) %>%
  ggplot(aes(p,cp)) +
  geom_point() +
  geom_line()

# Decidimos quedarnos con el de tres 
summary(pp_models)

## Subset selection object
## Call: regsubsets.formula(PE ~ ., data = power_plant, nvmax = 2500)
## 4 Variables  (and intercept)
##    Forced in Forced out
## AT     FALSE      FALSE
## V      FALSE      FALSE
## AP     FALSE      FALSE
## RH     FALSE      FALSE
## 1 subsets of each size up to 4
## Selection Algorithm: exhaustive
##          AT  V   AP  RH 
## 1  ( 1 ) "*" " " " " " "
## 2  ( 1 ) "*" " " " " "*"
## 3  ( 1 ) "*" "*" " " "*"
## 4  ( 1 ) "*" "*" "*" "*"

pp_model_cp <- lm(PE~AT+V+RH, data=power_plant)

anova(lm(PE~., data=power_plant))

## Analysis of Variance Table
## 
## Response: PE
##             Df  Sum Sq Mean Sq   F value    Pr(>F)    
## AT           1 2505095 2505095 120563.32 < 2.2e-16 ***
## V            1   46766   46766   2250.72 < 2.2e-16 ***
## AP           1    6259    6259    301.23 < 2.2e-16 ***
## RH           1   29875   29875   1437.81 < 2.2e-16 ***
## Residuals 9563  198702      21                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pp_model_anova <- lm(PE~., data=power_plant)

glance(pp_model_cp)

## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic p.value    df  logLik    AIC
## *     <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>   <dbl>  <dbl>
## 1     0.928         0.928  4.57    41321.       0     4 -28110. 56229.
## # ... with 3 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>

glance(pp_model_anova)

## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic p.value    df  logLik    AIC
## *     <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>   <dbl>  <dbl>
## 1     0.929         0.929  4.56    31138.       0     5 -28088. 56188.
## # ... with 3 more variables: BIC <dbl>, deviance <dbl>, df.residual <int>

anova(pp_model_cp, pp_model_anova)

## Analysis of Variance Table
## 
## Model 1: PE ~ AT + V + RH
## Model 2: PE ~ AT + V + AP + RH
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1   9564 199598                                  
## 2   9563 198702  1    895.28 43.087 5.507e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ambos modelos tienen $R^2$ y AIC similares, pero la prueba F declara mejor al modelo completo.

pp_model_data <- pp_model_anova %>%
  augment()

pp_model_data %>%
  ggplot(aes(x=.fitted, y=.resid)) +
  geom_point(alpha=0.5) +
  geom_smooth()

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Se alcanza a distinguir que en valores pequeños los residuales son mayores.

pp_model_data %>%
  ggplot(aes(1:9568, .cooksd)) +
  geom_point() +
  geom_line()

pp_model_data %>%
  ggplot(aes(1:9568, .hat)) +
  geom_point() +
  geom_line()

Ningún punto parece tener demasiado alto apalancamiento.

pp_model_data %>%
  ggplot(aes(sample=.std.resid)) +
  geom_qq_band(alpha=0.5, distribution='norm') +
    #geom_qq_band(alpha = 0.5, distribution='t', dparams=17) +
  stat_qq_line() +
  stat_qq_point() +
  labs(x = "Theoretical quantiles", y = "standarized residuals")

Pero los errores definitivamente no son normales, por lo que no podemos hacer inferencia sobre los coeficientes.

Haciendo la descomposición en componentes principales

pp_pc <- power_plant %>%
  select(-PE) %>%
  prcomp(scale=T) %>%
  print()

## Standard deviations (1, .., p=4):
## [1] 1.5615895 0.9533469 0.7416225 0.3202560
## 
## Rotation (n x k) = (4 x 4):
##           PC1         PC2       PC3        PC4
## AT  0.6148135 -0.04980763 0.1667322 -0.7692359
## V   0.5596829  0.09961674 0.5936655  0.5695549
## AP -0.4041859 -0.64834308 0.6287465 -0.1447856
## RH -0.3813045  0.75315799 0.4737504 -0.2508395

summary(pp_pc)

## Importance of components:
##                           PC1    PC2    PC3     PC4
## Standard deviation     1.5616 0.9533 0.7416 0.32026
## Proportion of Variance 0.6096 0.2272 0.1375 0.02564
## Cumulative Proportion  0.6096 0.8369 0.9744 1.00000

Aunque la primera componente principal es difícil de interpretar, la segunda tiene una interpretación que coincide con lo dicho en la documentación de los datos: separa vacumm —que afecta a las turbinas de vapor— , de las tres variables ambientales, que afectan a las de gas.

Aunque se necesitan tres de las cuatro para explicar la variabilidad de los predictores, vamos a ver cuántas se necesitan para explicar la de la respuesta.

pcr_model <- pcr(PE~., data=power_plant)

summary(pcr_model)

## Data:    X dimension: 9568 4 
##  Y dimension: 9568 1
## Fit method: svdpc
## Number of components considered: 4
## TRAINING: % variance explained
##     1 comps  2 comps  3 comps  4 comps
## X     62.99    92.27    98.41   100.00
## PE    62.17    84.83    85.79    92.87

Con la variabilidad de la respuesta nos va peor.

predplot(pcr_model)

validationplot(pcr_model, val.type = 'MSEP')

Sin embargo, en términos del error medio es parece que quedarnos con dos componentes es lo mejor. Hacemos la prueba formal

predictores_cp <- pp_pc$x %>%
  as_data_frame() %>%
  mutate(PE=power_plant$PE) %>%
  select(PE, everything()) %>% 
  print()

## # A tibble: 9,568 x 5
##       PE    PC1     PC2     PC3     PC4
##    <dbl>  <dbl>   <dbl>   <dbl>   <dbl>
##  1  463. -1.67  -1.25    0.449  -0.339 
##  2  444.  0.747 -1.44    0.784  -0.104 
##  3  489. -2.27   1.07   -0.527   0.536 
##  4  446.  0.351  0.517  -0.0437  0.0267
##  5  474. -1.80   1.57   -0.653  -0.144 
##  6  444.  1.22  -0.642  -0.193  -0.178 
##  7  467. -0.868 -0.0394 -0.424  -0.127 
##  8  478. -1.48  -1.00   -0.279   0.595 
##  9  476. -0.566 -2.62   -0.685   0.443 
## 10  478. -1.19  -0.370  -0.564   0.334 
## # ... with 9,558 more rows

names(predictores_cp) <- c('PE', 'PC1', 'PC2', 'PC3', 'PC4')

pcr_model_full <- lm(PE~.,data=predictores_cp)
  
summary(pcr_model_full)

## 
## Call:
## lm(formula = PE ~ ., data = predictores_cp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -43.435  -3.166  -0.118   3.201  17.778 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 454.36501    0.04660 9750.14   <2e-16 ***
## PC1          -9.99355    0.02984 -334.87   <2e-16 ***
## PC2          -1.53914    0.04888  -31.49   <2e-16 ***
## PC3          -5.08334    0.06284  -80.89   <2e-16 ***
## PC4          10.16892    0.14552   69.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.558 on 9563 degrees of freedom
## Multiple R-squared:  0.9287, Adjusted R-squared:  0.9287 
## F-statistic: 3.114e+04 on 4 and 9563 DF,  p-value: < 2.2e-16

pcr_model_full %>%
  augment() %>%
  ggplot(aes(sample=.std.resid)) +
  geom_qq_band(alpha=0.5, distribution='norm') +
  stat_qq_line() +
  stat_qq_point() +
  labs(x = "Theoretical quantiles", y = "standarized residuals")

Lamentablemente, no podemos hacer inferencia sobre los coeficientes porque el error no es normal.